Anaerobic CUE

Roey Angel 2021-03-02

Filter out rare taxa and those not classified as bacteria or archaea

Setting general parameters:

set.seed(2021)
prev_thresh <- 0.1 # remove ASVs with prevalence less than that
samples_prep_path <- "./"
data_path <- "./DADA2_pseudo/"
# Metadata_table <- "./AnCUE_Metadata_decontam.csv"
# Seq_table <- "DADA2.seqtab_nochim_decontam.tsv"
# Seq_file <- "DADA2.Seqs_decontam.fa"
Ps_file <- "Ps_obj_decontam.Rds"
Seq_file <- "DADA2.Seqs_decontam.fa"

Read phyloseq object

Also remove controls

readRDS(paste0(data_path, Ps_file)) ->
  Ps_obj

Ps_obj %>% 
  subset_samples(., Hours == 0) ->
  Ps_obj_t0

Ps_obj  %>% 
  subset_samples(., Hours != 0 & !is.na(Site)) -> 
  Ps_obj_SIP 

Combine repeated runs

sample_data(Ps_obj_SIP) %<>% 
  as("data.frame") %>% 
  rownames_to_column() %>% 
  mutate(Identifier = paste(Site, Oxygen, Glucose, Label..13C., Hours, Fraction.no., 
                            sep = "_")) %>% 
  column_to_rownames("rowname") 
  
phyloseq_merge_samples(Ps_obj_SIP, "Identifier") %>% 
  filter_taxa(., function(x) sum(x) > 0, TRUE) ->
  Ps_obj_SIP_merged 

# Compute lib sizes
sample_data(Ps_obj_SIP_merged)$Lib.size <- rowSums(otu_table(Ps_obj_SIP_merged))

Exploring dataset features

First let us look at the count data distribution

plot_lib_size(Ps_obj_SIP_merged, x = "Fraction.no.", fill = "Oxygen", facet1 = "Site", facet2 = "Hours")

I will test now the effect of library size and all other experimental factors on the community composition and also plot

(mod1 <- adonis(vegdist(otu_table(Ps_obj_SIP_merged), method = "bray") ~ Site * Oxygen * Hours + Lib.size,
  data = as(sample_data(Ps_obj_SIP_merged), "data.frame"),
  permutations = 999
))

## 
## Call:
## adonis(formula = vegdist(otu_table(Ps_obj_SIP_merged), method = "bray") ~      Site * Oxygen * Hours + Lib.size, data = as(sample_data(Ps_obj_SIP_merged),      "data.frame"), permutations = 999) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                    Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Site                1    22.555 22.5553 314.924 0.39126  0.001 ***
## Oxygen              1     6.377  6.3767  89.033 0.11061  0.001 ***
## Hours               1     0.333  0.3332   4.652 0.00578  0.002 ** 
## Lib.size            1     4.943  4.9430  69.016 0.08575  0.001 ***
## Site:Oxygen         1     1.483  1.4830  20.706 0.02573  0.001 ***
## Site:Hours          1     0.471  0.4707   6.572 0.00817  0.002 ** 
## Oxygen:Hours        1     0.281  0.2809   3.922 0.00487  0.005 ** 
## Site:Oxygen:Hours   1     0.363  0.3631   5.070 0.00630  0.001 ***
## Residuals         291    20.842  0.0716         0.36154           
## Total             299    57.648                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_lib_dist(Ps_obj_SIP_merged)

plot_read_dist(Ps_obj_SIP_merged)

plot_mean_SD(Ps_obj_SIP_merged)

Modelling library size shows a significant effect of read depth on the community structure, but explaining only 9% of the variance. The reads histogram shows as expected a highly sparse and skewed sequence matrix. The mean vs SD also shows as expected large dependency of SD on the mean reads of a sequence across all samples.

Taxa-based filtering

Now let us look at the taxonomic distribution

table(tax_table(Ps_obj_SIP_merged)[, "Kingdom"], exclude = NULL)

## 
##  Archaea Bacteria 
##      185    14280

table(tax_table(Ps_obj_SIP_merged)[, "Class"], exclude = NULL)

## 
##          Abditibacteria          Acidimicrobiia          Acidobacteriae 
##                      21                     679                    1199 
##          Actinobacteria                     AD3     Alphaproteobacteria 
##                     784                      27                    2031 
##            Anaerolineae           Armatimonadia                Babeliae 
##                       3                      84                     680 
##                 Bacilli             Bacteroidia                  BD7-11 
##                     459                     503                      16 
##         Bdellovibrionia          Blastocatellia         Campylobacteria 
##                     110                      20                       5 
##              Chlamydiae            Chloroflexia        Chthonomonadetes 
##                     583                       4                      63 
##              Clostridia          Coriobacteriia          Cyanobacteriia 
##                     276                       4                     178 
##         Dehalococcoidia              Deinococci      Desulfitobacteriia 
##                       1                      10                      12 
##           Desulfobulbia        Desulfovibrionia           Elusimicrobia 
##                       1                       2                      14 
##            Endomicrobia           Fibrobacteria          Fimbriimonadia 
##                       2                       5                      58 
##           Fusobacteriia     Gammaproteobacteria        Gemmatimonadetes 
##                      21                    2250                      25 
##             Gitt-GS-136         Gracilibacteria            Halobacteria 
##                       1                       1                       1 
##              Holophagae         Hydrogenedentia            JG30-KF-CM66 
##                       2                       1                      25 
##            Kapabacteria                  KD4-96         Ktedonobacteria 
##                      10                       1                      82 
##           Lentisphaeria             Lineage IIa           Longimicrobia 
##                       2                       5                       2 
##         Methanobacteria            Micrarchaeia              Myxococcia 
##                       1                       4                     111 
##           Negativicutes         Nitrososphaeria             Oligoflexia 
##                      13                     129                     106 
##                   OM190             Omnitrophia           Parcubacteria 
##                       1                      11                      11 
##           Phycisphaerae          Planctomycetes               Polyangia 
##                     154                    1703                     146 
##            Rhodothermia S0134 terrestrial group         Saccharimonadia 
##                       2                       1                      49 
##       Sericytochromatia                  SJA-28            Spirochaetia 
##                       1                       1                       1 
##             Subgroup 22              Subgroup 5       Syntrophobacteria 
##                       1                      12                       1 
##     Thermoanaerobaculia         Thermoleophilia          Thermoplasmata 
##                       2                     223                      50 
##                    TK10            Unclassified               vadinHA49 
##                      10                     615                      11 
##        Vampirivibrionia        Verrucomicrobiae        Vicinamibacteria 
##                     206                     599                      27

# table(tax_table(Ps_obj)[, "Family"], exclude = NULL)

Accordingly, we will remove some taxa which are obvious artefacts or those which aren’t bacteria or archaea

kingdoms2remove <- c("", "Eukaryota", "Unclassified")
orders2remove <- c("Chloroplast")
families2remove <- c("Mitochondria")

Ps_obj_kingdoms <- tax_glom(Ps_obj_SIP_merged, "Kingdom")
Ps_obj_orders <- tax_glom(Ps_obj_SIP_merged, "Order")
Ps_obj_families <- tax_glom(Ps_obj_SIP_merged, "Family")

Ps_obj_SIP_merged_filt <- subset_taxa(Ps_obj_SIP_merged, !is.na(Phylum) &
                        !Kingdom %in% kingdoms2remove &
                      !Order %in% orders2remove &
                      !Family %in% families2remove)

Summary_pruned <- tibble(
  Level = c("Kingdom", "Order", "Family"),
  ASVs.removed = c(
    table(tax_table(Ps_obj)[, "Kingdom"], exclude = NULL) %>% as.data.frame() %>% .[.$Var1 == "" | .$Var1 == "Eukaryota" | .$Var1 == "Unclassified", 2] %>% sum(),
    table(tax_table(Ps_obj)[, "Order"], exclude = NULL) %>% as.data.frame() %>% .[.$Var1 == "Chloroplast", 2] %>% sum(),
    table(tax_table(Ps_obj)[, "Family"], exclude = NULL) %>% as.data.frame() %>% .[.$Var1 == "Mitochondria", 2] %>% sum()
                     ),
  Seqs.removed = c(
    psmelt(Ps_obj_kingdoms) %>%
      group_by(Kingdom) %>%
      filter(Kingdom == "" |
               Kingdom == "Eukaryota" | Kingdom == "Unclassified") %>%
      summarise(sum = sum(Abundance)) %>% .$sum %>% sum(),
    psmelt(Ps_obj_orders) %>%
      group_by(Order) %>%
      filter(Order == orders2remove) %>%
      summarise(sum = sum(Abundance)) %>% .$sum %>% sum(),
    psmelt(Ps_obj_families) %>%
      group_by(Family) %>%
      filter(Family == families2remove) %>%
      summarise(sum = sum(Abundance)) %>% .$sum %>% sum()
    )
  )

Summary_pruned %>% 
  kable(., digits = c(0, 1, 0)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F)

Removed 0.0335% of the sequences.

Now let’s explore the prevalence of different taxa in the database. Prevalence is the number of samples in which a taxa appears at least once. So “Mean prevalence” refers to in how many samples does a sequence belonging to the phylum appears on average, and “Sum prevalence” is the sum of all samples where any sequence from the taxon appears.

prevdf <- apply(X = otu_table(Ps_obj_SIP_merged_filt),
                 MARGIN = ifelse(taxa_are_rows(Ps_obj_SIP_merged_filt), yes = 1, no = 2),
                 FUN = function(x){sum(x > 0)})
# Add taxonomy and total read counts to this data.frame
prevdf <- data.frame(Prevalence = prevdf,
                      TotalAbundance = taxa_sums(Ps_obj_SIP_merged_filt),
                      tax_table(Ps_obj_SIP_merged_filt))

prevdf %>%
  group_by(Phylum) %>%
  summarise(`Mean prevalence` = mean(Prevalence),
            `Sum prevalence` = sum(Prevalence)) ->
  Prevalence_phylum_summary

Prevalence_phylum_summary %>% 
  kable(., digits = c(0, 1, 0)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F)

prevdf %>%
  group_by(Order) %>%
  summarise(`Mean prevalence` = mean(Prevalence),
            `Sum prevalence` = sum(Prevalence)) ->
  Prevalence_order_summary

Prevalence_order_summary %>% 
  kable(., digits = c(0, 1, 0)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F)

Based on that I’ll remove all orders with a sum prevalence of under 5% (15) of all samples

Prevalence_order_summary %>% 
  filter(`Sum prevalence` < (0.05 * nsamples(Ps_obj_SIP_merged_filt))) %>% 
  dplyr::select(Order) %>% 
  map(as.character) %>% 
  unlist() ->
  filterOrder

Ps_obj_SIP_merged_filt2 <- subset_taxa(Ps_obj_SIP_merged_filt, !Order %in% filterOrder)

sample_data(Ps_obj_SIP_merged_filt2)$Lib.size <- rowSums(otu_table(Ps_obj_SIP_merged_filt2))
print(Ps_obj_SIP_merged_filt)

## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 14199 taxa and 300 samples ]:
## sample_data() Sample Data:        [ 300 samples by 29 sample variables ]:
## tax_table()   Taxonomy Table:     [ 14199 taxa by 6 taxonomic ranks ]:
## taxa are columns

print(Ps_obj_SIP_merged_filt2)

## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 14099 taxa and 300 samples ]:
## sample_data() Sample Data:        [ 300 samples by 29 sample variables ]:
## tax_table()   Taxonomy Table:     [ 14099 taxa by 6 taxonomic ranks ]:
## taxa are columns

This removed 100 or 1% of the ESVs, and 0.046% of the reads.

Plot general prevalence features of the phyla

# Subset to the remaining phyla
prevdf_phylum_filt <- subset(prevdf, 
                             Phylum %in% get_taxa_unique(Ps_obj_SIP_merged_filt2, 
                                                         "Phylum"))
ggplot(prevdf_phylum_filt,
       aes(TotalAbundance, 
           Prevalence / nsamples(Ps_obj_SIP_merged_filt2), 
           color = Phylum)) +
  # Include a guess for parameter
  geom_hline(yintercept = 0.05,
             alpha = 0.5,
             linetype = 2) + 
  geom_point2(size = 2, alpha = 0.7) +
  scale_x_log10() +  
  xlab("Total Abundance") + 
  ylab("Prevalence [Frac. Samples]") +
  facet_wrap( ~ Phylum) + 
  theme(legend.position = "none")

Plot general prevalence features of the top 20 orders

# Subset to the remaining phyla
prevdf_order_filt <- subset(prevdf, 
                            Order %in% get_taxa_unique(Ps_obj_SIP_merged_filt2, "Order"))

# grab the top 30 most abundant orders
prevdf_order_filt %>% 
  group_by(Order) %>%
  summarise(Combined.abundance = sum(TotalAbundance)) %>% 
  arrange(desc(Combined.abundance)) %>% 
  .[1:30, "Order"]  ->
  Orders2plot

prevdf_order_filt2 <- subset(prevdf,
                             Order %in% Orders2plot$Order)

ggplot(prevdf_order_filt2,
       aes(TotalAbundance,
           Prevalence / nsamples(Ps_obj_SIP_merged_filt2), 
           color = Order)) +
  # Include a guess for parameter
  geom_hline(yintercept = 0.05,
             alpha = 0.5,
             linetype = 2) + 
  geom_point2(size = 2, alpha = 0.7) +
  scale_x_log10() +  
  xlab("Total Abundance") +
  ylab("Prevalence [Frac. Samples]") +
  facet_wrap( ~ Order) + 
  theme(legend.position = "none")

Unsupervised filtering by prevalence

I’ll remove all sequences which appear in less than 5% of the samples

# Define prevalence threshold as 5% of total samples
prevalenceThreshold <- prev_thresh * nsamples(Ps_obj_SIP_merged_filt2)
prevalenceThreshold

## [1] 30

# Execute prevalence filter, using `prune_taxa()` function
keepTaxa <-
  row.names(prevdf_phylum_filt)[(prevdf_phylum_filt$Prevalence >= prevalenceThreshold)]
Ps_obj_SIP_merged_filt3 <- prune_taxa(keepTaxa, Ps_obj_SIP_merged_filt2)

sample_data(Ps_obj_SIP_merged_filt3)$Lib.size <- rowSums(otu_table(Ps_obj_SIP_merged_filt3))
print(Ps_obj_SIP_merged_filt2)

## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 14099 taxa and 300 samples ]:
## sample_data() Sample Data:        [ 300 samples by 29 sample variables ]:
## tax_table()   Taxonomy Table:     [ 14099 taxa by 6 taxonomic ranks ]:
## taxa are columns

print(Ps_obj_SIP_merged_filt3)

## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 3693 taxa and 300 samples ]:
## sample_data() Sample Data:        [ 300 samples by 29 sample variables ]:
## tax_table()   Taxonomy Table:     [ 3693 taxa by 6 taxonomic ranks ]:
## taxa are columns

This removed 10406 or 74% of the ESVs!

However all these removed ESVs accounted for only:

prevdf_phylum_filt %>% 
  arrange(., Prevalence) %>% 
  group_by(Prevalence > prevalenceThreshold) %>% 
  summarise(Abundance = sum(TotalAbundance)) %>%
  mutate(`Rel. Ab.` = percent(Abundance / sum(Abundance), accuracy = 0.01)) %>% 
  kable(., digits = c(0, 1, 0)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F)

So it’s fine to remove them.

Test again the effect of library size and all other experimental factors on the community composition after filtering

(mod1 <- adonis(vegdist(otu_table(Ps_obj_SIP_merged_filt3), method = "bray") ~ Site * Oxygen * Hours + Lib.size,
  data = as(sample_data(Ps_obj_SIP_merged_filt3), "data.frame"),
  permutations = 999
))

## 
## Call:
## adonis(formula = vegdist(otu_table(Ps_obj_SIP_merged_filt3),      method = "bray") ~ Site * Oxygen * Hours + Lib.size, data = as(sample_data(Ps_obj_SIP_merged_filt3),      "data.frame"), permutations = 999) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                    Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Site                1    22.638 22.6379  342.15 0.40294  0.001 ***
## Oxygen              1     6.358  6.3575   96.09 0.11316  0.001 ***
## Hours               1     0.322  0.3223    4.87 0.00574  0.005 ** 
## Lib.size            1     5.069  5.0686   76.61 0.09022  0.001 ***
## Site:Oxygen         1     1.475  1.4747   22.29 0.02625  0.001 ***
## Site:Hours          1     0.458  0.4576    6.92 0.00814  0.001 ***
## Oxygen:Hours        1     0.264  0.2643    3.99 0.00470  0.007 ** 
## Site:Oxygen:Hours   1     0.345  0.3448    5.21 0.00614  0.002 ** 
## Residuals         291    19.254  0.0662         0.34270           
## Total             299    56.181                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_lib_dist(Ps_obj_SIP_merged_filt3)

plot_read_dist(Ps_obj_SIP_merged_filt3)

plot_mean_SD(Ps_obj_SIP_merged_filt3)

Save filtered phyloseq object

saveRDS(Ps_obj_SIP_merged_filt3, file = paste0(data_path, str_remove(Ps_file, ".Rds"), "_filt3.Rds"))
readDNAStringSet(paste0(data_path, Seq_file)) %>% 
  .[taxa_names(Ps_obj_SIP_merged_filt3)] %>%  
  writeXStringSet(., filepath = paste0(data_path, str_remove(Seq_file, ".fa*"), "_filtered.fa"), format = "fasta", width = 1000)

sessioninfo::session_info() %>%
  details::details(
    summary = 'Current session info',
    open    = TRUE
  )

─ Session info ─────────────────────────────────────────────────────────────────────────
 setting  value                       
 version  R version 4.0.3 (2020-10-10)
 os       Ubuntu 18.04.5 LTS          
 system   x86_64, linux-gnu           
 ui       X11                         
 language (EN)                        
 collate  en_US.UTF-8                 
 ctype    en_US.UTF-8                 
 tz       Europe/Prague               
 date     2021-03-02                  

─ Packages ─────────────────────────────────────────────────────────────────────────────
 package        * version    date       lib source                           
 ade4             1.7-16     2020-10-28 [1] CRAN (R 4.0.2)                   
 affy             1.66.0     2020-04-27 [1] Bioconductor                     
 affyio           1.58.0     2020-04-27 [1] Bioconductor                     
 ape              5.4-1      2020-08-13 [1] CRAN (R 4.0.2)                   
 assertthat       0.2.1      2019-03-21 [1] CRAN (R 4.0.2)                   
 backports        1.2.1      2020-12-09 [1] CRAN (R 4.0.2)                   
 bayestestR       0.8.2      2021-01-26 [1] CRAN (R 4.0.3)                   
 Biobase        * 2.48.0     2020-04-27 [1] Bioconductor                     
 BiocGenerics   * 0.34.0     2020-04-27 [1] Bioconductor                     
 BiocManager      1.30.10    2019-11-16 [1] CRAN (R 4.0.2)                   
 biomformat       1.16.0     2020-04-27 [1] Bioconductor                     
 Biostrings     * 2.56.0     2020-04-27 [1] Bioconductor                     
 broom            0.7.5      2021-02-19 [1] CRAN (R 4.0.3)                   
 cellranger       1.1.0      2016-07-27 [1] CRAN (R 4.0.2)                   
 cli              2.3.1      2021-02-23 [1] CRAN (R 4.0.3)                   
 clipr            0.7.1      2020-10-08 [1] CRAN (R 4.0.2)                   
 cluster          2.1.1      2021-02-14 [1] CRAN (R 4.0.3)                   
 coda             0.19-4     2020-09-30 [1] CRAN (R 4.0.2)                   
 codetools        0.2-18     2020-11-04 [1] CRAN (R 4.0.2)                   
 colorspace       2.0-0      2020-11-11 [1] CRAN (R 4.0.2)                   
 crayon           1.4.1      2021-02-08 [1] CRAN (R 4.0.3)                   
 data.table       1.14.0     2021-02-21 [1] CRAN (R 4.0.3)                   
 DBI              1.1.1      2021-01-15 [1] CRAN (R 4.0.3)                   
 dbplyr           2.1.0      2021-02-03 [1] CRAN (R 4.0.3)                   
 desc             1.2.0      2018-05-01 [1] CRAN (R 4.0.2)                   
 details          0.2.1      2020-01-12 [1] CRAN (R 4.0.2)                   
 digest           0.6.27     2020-10-24 [1] CRAN (R 4.0.2)                   
 dplyr          * 1.0.4      2021-02-02 [1] CRAN (R 4.0.3)                   
 effectsize       0.4.3      2021-01-18 [1] CRAN (R 4.0.3)                   
 ellipsis         0.3.1      2020-05-15 [1] CRAN (R 4.0.2)                   
 emmeans          1.5.4      2021-02-03 [1] CRAN (R 4.0.3)                   
 estimability     1.3        2018-02-11 [1] CRAN (R 4.0.2)                   
 evaluate         0.14       2019-05-28 [1] CRAN (R 4.0.2)                   
 extrafont      * 0.17       2014-12-08 [1] CRAN (R 4.0.2)                   
 extrafontdb      1.0        2012-06-11 [1] CRAN (R 4.0.2)                   
 fansi            0.4.2      2021-01-15 [1] CRAN (R 4.0.3)                   
 farver           2.1.0      2021-02-28 [1] CRAN (R 4.0.3)                   
 forcats        * 0.5.1      2021-01-27 [1] CRAN (R 4.0.3)                   
 foreach          1.5.1      2020-10-15 [1] CRAN (R 4.0.2)                   
 fs               1.5.0      2020-07-31 [1] CRAN (R 4.0.2)                   
 generics         0.1.0      2020-10-31 [1] CRAN (R 4.0.2)                   
 ggplot2        * 3.3.3      2020-12-30 [1] CRAN (R 4.0.2)                   
 ggridges         0.5.3      2021-01-08 [1] CRAN (R 4.0.2)                   
 glue             1.4.2      2020-08-27 [1] CRAN (R 4.0.2)                   
 gtable           0.3.0      2019-03-25 [1] CRAN (R 4.0.2)                   
 haven            2.3.1      2020-06-01 [1] CRAN (R 4.0.2)                   
 hexbin           1.28.2     2021-01-08 [1] CRAN (R 4.0.2)                   
 highr            0.8        2019-03-20 [1] CRAN (R 4.0.2)                   
 hms              1.0.0      2021-01-13 [1] CRAN (R 4.0.3)                   
 htmltools        0.5.1.1    2021-01-22 [1] CRAN (R 4.0.3)                   
 httr             1.4.2      2020-07-20 [1] CRAN (R 4.0.2)                   
 igraph           1.2.6      2020-10-06 [1] CRAN (R 4.0.2)                   
 insight          0.13.1     2021-02-22 [1] CRAN (R 4.0.3)                   
 IRanges        * 2.22.2     2020-05-21 [1] Bioconductor                     
 iterators        1.0.13     2020-10-15 [1] CRAN (R 4.0.2)                   
 jsonlite         1.7.2      2020-12-09 [1] CRAN (R 4.0.2)                   
 kableExtra     * 1.3.4      2021-02-20 [1] CRAN (R 4.0.3)                   
 knitr            1.31       2021-01-27 [1] CRAN (R 4.0.3)                   
 labeling         0.4.2      2020-10-20 [1] CRAN (R 4.0.2)                   
 lattice        * 0.20-41    2020-04-02 [1] CRAN (R 4.0.2)                   
 lifecycle        1.0.0      2021-02-15 [1] CRAN (R 4.0.3)                   
 limma            3.44.3     2020-06-12 [1] Bioconductor                     
 lubridate        1.7.10     2021-02-26 [1] CRAN (R 4.0.3)                   
 magrittr       * 2.0.1      2020-11-17 [1] CRAN (R 4.0.2)                   
 MASS             7.3-53.1   2021-02-12 [1] CRAN (R 4.0.3)                   
 Matrix           1.3-2      2021-01-06 [1] CRAN (R 4.0.2)                   
 mgcv             1.8-34     2021-02-16 [1] CRAN (R 4.0.3)                   
 modelr           0.1.8      2020-05-19 [1] CRAN (R 4.0.2)                   
 multcomp         1.4-16     2021-02-08 [1] CRAN (R 4.0.3)                   
 multtest         2.44.0     2020-04-27 [1] Bioconductor                     
 munsell          0.5.0      2018-06-12 [1] CRAN (R 4.0.2)                   
 mvtnorm          1.1-1      2020-06-09 [1] CRAN (R 4.0.2)                   
 nlme             3.1-152    2021-02-04 [1] CRAN (R 4.0.3)                   
 parameters       0.12.0     2021-02-21 [1] CRAN (R 4.0.3)                   
 permute        * 0.9-5      2019-03-12 [1] CRAN (R 4.0.2)                   
 phyloseq       * 1.32.0     2020-04-27 [1] Bioconductor                     
 pillar           1.5.0      2021-02-22 [1] CRAN (R 4.0.3)                   
 pkgconfig        2.0.3      2019-09-22 [1] CRAN (R 4.0.2)                   
 plyr             1.8.6      2020-03-03 [1] CRAN (R 4.0.2)                   
 png              0.1-7      2013-12-03 [1] CRAN (R 4.0.2)                   
 preprocessCore   1.50.0     2020-04-27 [1] Bioconductor                     
 prettyunits      1.1.1      2020-01-24 [1] CRAN (R 4.0.2)                   
 progress         1.2.2      2019-05-16 [1] CRAN (R 4.0.2)                   
 purrr          * 0.3.4      2020-04-17 [1] CRAN (R 4.0.2)                   
 R6               2.5.0      2020-10-28 [1] CRAN (R 4.0.2)                   
 RColorBrewer     1.1-2      2014-12-07 [1] CRAN (R 4.0.2)                   
 Rcpp             1.0.6      2021-01-15 [1] CRAN (R 4.0.3)                   
 readr          * 1.4.0      2020-10-05 [1] CRAN (R 4.0.2)                   
 readxl           1.3.1      2019-03-13 [1] CRAN (R 4.0.2)                   
 reprex           1.0.0      2021-01-27 [1] CRAN (R 4.0.3)                   
 reshape2         1.4.4      2020-04-09 [1] CRAN (R 4.0.2)                   
 rhdf5            2.32.4     2020-10-05 [1] Bioconductor                     
 Rhdf5lib         1.10.1     2020-07-09 [1] Bioconductor                     
 rlang            0.4.10     2020-12-30 [1] CRAN (R 4.0.2)                   
 rmarkdown        2.7        2021-02-19 [1] CRAN (R 4.0.3)                   
 rprojroot        2.0.2      2020-11-15 [1] CRAN (R 4.0.2)                   
 rstudioapi       0.13       2020-11-12 [1] CRAN (R 4.0.2)                   
 Rttf2pt1         1.3.8      2020-01-10 [1] CRAN (R 4.0.2)                   
 rvest            0.3.6      2020-07-25 [1] CRAN (R 4.0.2)                   
 S4Vectors      * 0.26.1     2020-05-16 [1] Bioconductor                     
 sandwich         3.0-0      2020-10-02 [1] CRAN (R 4.0.2)                   
 scales         * 1.1.1      2020-05-11 [1] CRAN (R 4.0.2)                   
 see            * 0.6.2      2021-02-04 [1] CRAN (R 4.0.3)                   
 sessioninfo      1.1.1      2018-11-05 [1] CRAN (R 4.0.2)                   
 speedyseq      * 0.5.3.9001 2020-10-27 [1] Github (mikemc/speedyseq@8daed32)
 stringi          1.5.3      2020-09-09 [1] CRAN (R 4.0.2)                   
 stringr        * 1.4.0      2019-02-10 [1] CRAN (R 4.0.2)                   
 survival         3.2-7      2020-09-28 [1] CRAN (R 4.0.2)                   
 svglite        * 2.0.0      2021-02-20 [1] CRAN (R 4.0.3)                   
 systemfonts      1.0.1      2021-02-09 [1] CRAN (R 4.0.3)                   
 TH.data          1.0-10     2019-01-21 [1] CRAN (R 4.0.2)                   
 tibble         * 3.1.0      2021-02-25 [1] CRAN (R 4.0.3)                   
 tidyr          * 1.1.2      2020-08-27 [1] CRAN (R 4.0.2)                   
 tidyselect       1.1.0      2020-05-11 [1] CRAN (R 4.0.2)                   
 tidyverse      * 1.3.0      2019-11-21 [1] CRAN (R 4.0.2)                   
 utf8             1.1.4      2018-05-24 [1] CRAN (R 4.0.2)                   
 vctrs            0.3.6      2020-12-17 [1] CRAN (R 4.0.2)                   
 vegan          * 2.5-7      2020-11-28 [1] CRAN (R 4.0.2)                   
 viridisLite      0.3.0      2018-02-01 [1] CRAN (R 4.0.2)                   
 vsn            * 3.56.0     2020-04-27 [1] Bioconductor                     
 webshot          0.5.2      2019-11-22 [1] CRAN (R 4.0.2)                   
 withr            2.4.1      2021-01-26 [1] CRAN (R 4.0.3)                   
 xfun             0.21       2021-02-10 [1] CRAN (R 4.0.3)                   
 xml2             1.3.2      2020-04-23 [1] CRAN (R 4.0.2)                   
 xtable           1.8-4      2019-04-21 [1] CRAN (R 4.0.2)                   
 XVector        * 0.28.0     2020-04-27 [1] Bioconductor                     
 yaml             2.2.1      2020-02-01 [1] CRAN (R 4.0.2)                   
 zlibbioc         1.34.0     2020-04-27 [1] Bioconductor                     
 zoo              1.8-8      2020-05-02 [1] CRAN (R 4.0.2)                   

[1] /home/angel/R/library
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library

References